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Abstract 

We present a new three-dimensional Monte-Carlo code MUSIC (MUon Simula- 
tion Code) for muon propagation through the rock. All processes of muon interaction 
with matter with high energy loss (including the knock-on electron production) are 
treated as stochastic processes. The angular deviation and lateral displacement of 
muons due to multiple scattering, as well as bremsstrahlung, pair production and 
inelastic scattering are taken into account. The code has been applied to obtain the 
energy distribution and angular and lateral deviations of single muons at different 
depths underground. The muon multiplicity distributions obtained with MUSIC 
and CORSIKA (Extensive Air Shower simulation code) are also presented. We dis- 
cuss the systematic uncertainties of the results due to different muon bremsstrahlung 
cross-sections. 



PACS number (s): 94.40T 

Keywords: High energy muons, Muons underground, Monte Carlo methods. Inter- 
action of particles with matter. 
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1 Introduction 



Muon transport through the rock plays an important role in 'underground physics', 
in particular in the underground cosmic-ray experiments. The muon 'depth-intensity 
curve' related to the muon energy spectrum at sea level and to the primary cosmic- 
ray specrum, muon multiplicity distribution and 'decoherence curve' related to the 
mass composition of primaries and to the characteristics of high-energy hadron- 
nucleus interactions are influenced by the muon propagation through the rock. The 
large underground detectors designed to search for rare events (such as proton decay, 
neutrino interactions etc.) are not free from the background of cosmic-ray muons 
and muon-produced secondary particles. The propagation of muons produced by 
high-energy astrophysical and atmosheric neutrinos in the rock should be also taken 
into account. These examples show the importance of an accurate three-dimensional 
simulations of the muon transport in matter. 

There are several three-dimensional codes for the simulation of muon propaga- 
tion through the rock (see, for example, [|l], ^ |[). They have been used to obtain the 
characteristics of muon flux and muon bundles deep underground and to describe the 
experimental data. However, the increase of the accuracy of the recent and future 
measurements requests an adequate increase of the accuracy of the calculations. The 
accuracy of the calculations is restricted by the uncertainties in the cross-sections 
and the assumptions made to simplify the computational procedure (due to the re- 
stricted CPU time). Recently, a new corrections to the muon bremsstrahlung and 
knock-on electron production cross-sections were calculated [Q. A new generation 
of powerful computers allows to avoid many simplifications in the calculation pro- 
cedure. 

In our code, named MUSIC (MUon Simulation Code) we have used the most 
recent and accurate cross-sections of the muon interactions with matter and we 
tried to avoid many simplifications used in the previous simulations. This three- 
dimensional code is based on an algorithm of a one-dimensional muon propagation 
code described in and used to fit the 'depth-intensity relation' measured by LVD 
1^, 0]. We have treated all processes of muon interaction with high energy transfer 
(bremsstrahlung, inelastic scattering, pair production and knock-on electron pro- 
duction) as the stochastic ones if the fraction of energy lost by muon is more than 
10~^. The mean muon path between two interactions with the energy loss more than 
10~^ is about 20 hg/cm^. The angular deviation of muons was taken into account 
not only in the process of multiple scattering (as it was usually done) but also in 
the processes of bremsstrahlung, inelastic scattering and pair production. We inves- 
tigated the effect due to the last three processes. 

We have performed three-dimensional simulations of the muon transport thro- 
ugh Gran Sasso and standard rocks. The results of such simulations for Gran Sasso 
rock, as well as the complete description of the code and related topics can be found 
in p. Here we present the results of the calculations for standard rock using two 
different bremsstrahlung cross-sections ^. In this way we investigated the effect 
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due to the uncertainties in the cross-sections. We have also apphed the code to cal- 
culate muon multiplicity distribution using CORSIKA code for the simulations 
of Extensive Air Showers. 



The main features of the MUSIC code are described in Section 2. The results of 
the simulation in the standard rock are presented in Section 3. We also discuss in 
Section 3 the dependence of our results on the cross-sections and assumptions used 
in the code. The results of the full Monte-Carlo simulation of muon component 
using CORSIKA and MUSIC can be found in Section 4. In Section 5 we present 
our conclusions. 



In the code we have treated energy losses due to ionization, bremsstrahlung, pair 
production and muon - nucleus inelastic scattering. 
The mean energy loss is usually expressed by 



where a is the energy loss due to ionization and (3 is the relative energy loss due to 
other processes. 

Ionization can be considered as a continuous process and it is well described by 
Bethe-Bloch formula. On the contrary, in other processes a particle can lose a big 
fraction of its energy in a single interaction and it is necessary to treat these kinds 
of energy losses in two separate parts, one continuous and the other stochastic. To 
improve the simulation procedure and to study the importance of the stochasticity of 
the knock-on electron production we have divided the energy loss due to this process 
into two parts: continuous loss with v < Vcut = 10~^ which we calculate using the 
Bethe-Bloch formula, and stochastic loss which we treat in the same way as other 
stochastic processes. Thus, the energy losses due to knock-on electron production, 
bremsstrahlung, pair production and inelastic scattering are expressed by: 



where A is the mass number of the material, N is the Avogadro number, v is the 
fraction of the transferred muon energy and ^ is the cross-section of the process. 

The differential probability that a muon loses a fraction v of its energy per g/ cm^ 
in a stochastic way is given by 



2 The code 



- < — >= a{E)+l3{E)E 






It is quite difficult to choose a suitable value for Vcut, because a little value of 
Vcut, that means an accurate simulation, results in an increase of CPU's time. We 
have chosen Vcut = 10^'^ that seems meet all requirements. On a DEC AlphaStation 
250 4/266 the simulation of 10 TeV muons propagated through 3 km.w.e. takes 0.02 
s/muon. 

We note that most of the other codes use higher value of v < Vcut and/or are 
considering the knock-on electron production as a continuous process. 

The aim of the code is to propagate a muon through a layer X of standard 
rock (we used A=22, Z=ll, p= 2.65 g/cm^). If we know the initial energy of the 
particle Eq , its starting point Xq, z/q, Zq in the rock and its direction 6*0, (po, we can 
simulate with an iterative procedure, step by step, muon interaction point according 
to P(|T^|) = ^e~l ^ IZ-^W, where "r^ is the vector which defines the interaction point 
in the space and A is the mean free path. Then the continuous and stochastic 
energy losses, as well as lateral and angular displacement due to multiple Coulomb 
scattering are evaluated. Moreover the code simulates the angular deviations caused 
by the stochastic processes. 

Finally we check if the muon has crossed the layer X, and in this case we store 
its final energy Ej , its direction 6f, (pf and the coordinates of its arrival point 

To describe pair production we have used the cross-section given by Kokoulin- 
Petrukhin [ill]. The code samples separately muon energy loss and angular deflec- 



tion. The second one is sampled following the parametrization proposed in |T2| (see 
Appendix A). 

Two different muon bremsstrahlung cross-sections, one given by Bezrukov and 
Bugaev and the other recently proposed by Kelner, Kokoulin and Petrukhin 
1^ are available in MUSIC. These two calculations differ mainly by the account 
of bremsstrahlung on atomic electrons, when the photon is emitted by the target 
electron. 

The code samples the fraction of muon energy transferred to the photon and 
separately the muon scattering angle, following also in this case the parametrisation 
proposed by [0 (see Appendix A). 

We carried out additional test of the simulation of muon deflection angle due to 
pair production and bremsstrahlung process using a simplified procedure based on 



the GEANT parametrizations |T^, |T^. We concluded that this approach leads to a 
little underestimation of the muon deflection due to these processes with respect to 
Van Ginneken's parametrizations. 

Among stochastic processes, the angular deflection due to inelastic scattering on 



nuclei is expected to be dominant . In this case to calculate the muon energy loss 
and angular deflection we have used double-differential cross section jf^Q2 obtained 
by Bezrukov and Bugaev [|l6l, where Q is the four-momentum transferred. We have 
sampled the muon scattering angle according to the aforementioned cross-section 
using the relation = 2{E^QE^f — P^oP^fCosO^ — m?^). Since the muon inelastic 
scattering is expected to be the most important for angular deviation of muons (with 
the exception of multiple scattering), this process has been treated in a more correct 
way (using the double-differential cross-section) in comparison with other stochastic 



4 



processes. However, the cross-section of this process proposed in [0 is vahd only at 
low (the accuracy is not worse than 11% for <2> GeV^ and energy transfers 
from 1 to 10^ GeV 0). 



Some authors |T^, [T^, |T8| , pointed out the complexity of the problem to sample 
simultaneously the emission angle of the photon, the deflection angle of the muon and 
the fractional energy loss v. The accurate and simultaneous simulation of the angles 
and energies can be important for muons exceeding 1 TeV, where the deflection due 
to bremsstrahlung could be dominant [l^ . As it will be discussed in the next section. 



however, even if we used the same approach as in , we didn't find such behaviour, 
and the multiple scattering appears to be the dominant mechanism of deflection of 
muons passing through the large thickness of matter. A complete treatment of these 
aspects is currently beyond the aim of our code. 

We have treated multiple Coulomb scattering in the Gaussian approximation 
|2^. We have considered the projections of muon angular deviation, 6*^, and lateral 
displacement, /, on two orthogonal planes which include the initial direction of 
muon. The projections of angular deviation and lateral displacement of muon with 
respect to its initial direction at each step between two interactions with v > Vcut 
are distributed according to correlated Gaussians with mean values and dispersions 
calculated using the formulae: 



< >=< By >= 0, a', =< el >=< el >=^l% (4) 



< >=< >= =< ii >=< ii >= 1^ /; ^ (5) 

where e^s = 21.2 MeV, A = 26.48 g/cm^ is the radiation length in standard rock 
and p = 2.65 g/cm? is the standard rock density. As discussed in the next section, 
we have checked also if the use of a more accurate treatment of multiple scattering 
(e.g. Moliere theory) leads to different results. 



3 Results 

We have propagated 100000 muons through 3 km.w.e. and 1000000 through 10 km. 
w.e. with MUSIC using Kelner-Kokoulin-Petrukhin's cross-section for bremsstra- 
hlung. Initial muon energy has been sampled according to a primary energy spec- 
trum with index 7 = 3.7. In the case of 3 km.w.e. minimal muon energy Eq was 
fixed to 900 GeV, while for 10 km.w.e. Eq was 9 TeV, corresponding in each case 
to a survival probability of less or about 0.3%. 

Figures |l|, || and ^ show energy, lateral and angular distributions at 3 and 10 
km.w.e. The mean muon energy is 250 GeV at 3 km.w.e. and 367 GeV at 10 km.w.e. 
This result is in agreement with the simple consideration that the mean muon energy 
at a depth h increases with h. 

Mean angular deviation < a > is 0.56° at 3 km.w.e. and 0.45° at 10 km.w.e. 
Angular deviation is mainly caused by multiple Coulomb scattering. This process 
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produces a deviation that increases when energy decreases. This explains why at 3 
km.w.e., when < Eh > is lower, < a > is higher than at 10 km.w.e. 

Figures ^ and ^ show lateral and angular distributions obtained with MUSIC 
and PROPMU - a code for muon propagation in rock written by P. Lipari and T. 
Stanev 0. The codes are completely independent, but their results are practically 
the same. This fact confirms the correctness of both simulations and shows that the 
effects of angular deviations due to stochastic processes, with the approximations 
discussed in previous section, are really very little, 1-2 orders of magnitude less 
than the effects due to multiple scattering. In fact, as Lipari and Stanev treat only 
angular deviation due to multiple scattering, and in MUSIC all possible angular 
deviations are taken into account, the distributions obtained with these two codes 
do not show any appreciable difference. 

In Figure ^ we show the contribution of each process to the angular deviation 
of muons which pass 3 km.w.e. To do this plot we have repeated five times the 
same simulation leaving active one mechanism of deflection each time. For multiple 
scattering we have used either Gaussian treatment described in previous section or 
Moliere theory subroutines provided by GEANT [|19]. A more detailed treatment 
(Moliere theory) results in the similar distribution compared with that of Gaussian 
treatment, with a small increase of mean deflection angle. However, the original 
Moliere theory (see [|l^ and references therein) does not provide the lateral dis- 



placement, more important from the experimental point of view. As can be clearly 
seen from Figure ^ the multiple scattering dominates over other processes. This re- 
sult remains true for 10 km.w.e. and for a monochromatic beam of 10 TeV muons. 

To calculate survival probabilities, P{Eq, h), and muon energy distributions un- 
derground, P{Eq, h, E), we have propagated muons with different initial energy 
(100000 muons for each initial energy) through the standard rock. P{Eo, h, E) 
means the probability that muon with initial energy Eq will reach the depth h with 
energy E. Figure |^ shows the survival propabilities obtained with MUSIC for Eq 
from 1 to 10^ TeV. Figure § shows P{Eq, h) for h = 4 km.w.e. and h = 10 km.w.e., 
calculated with MUSIC and with PROPMU. The curves are very similar also in this 
case. Our values are on average (1-2)% lower: this difference can be due to the new 
cross-section for bremsstrahlung used in MUSIC. 

We have studied the effects of the cross-section difference and of the stochasticity 
of different processes on the characteristics of the muon flux underground. Using 
the sea-level vertical muon energy spectrum, proposed by T.Gaisser with the 
spectral index of primary spectrum 7 = 2.70, and the distributions P{EQ,h) and 
P{Eq, h, E) we have calculated the vertical muon intensities I{h) and mean muon 
energies < E{h) > at different depths underground: 

I{h)= rp{Eo,h)-^^-dE„ (6) 
JO aEn 



roo dI(E,h) ,rp 
Jo dE ' "-^ 



< m >= (7) 
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where is the vertical muon energy spectrum at sea level and '^^^J^^^ is the 

vertical muon spectrum at the depth h. The latter one is obtained using the formula: 

We have calculated the muon intensities and mean muon energies for two different 
muon bremsstrahlung cross-sections (from and 0) and for two values of Vcut 
(10^^ and 10~^). We have carried out additional tests considering the knock-on 
electron production and pair production as continuous processes. The results of 
the calculations are presented in Table |l] for I{h) and Table |^ for < E{h) >. The 
first column in both tables shows the depth in km.w.e. of standard rock. Other 
columns show the intensity in units cm~'^s~^sr~^ (Table |l]) and the mean muon 
energy in GeV (Table 2nd column present our basic values calculated with the 
muon bremsstrahlung cross-section from and Vcut = 10~^. The results obtained 
with bremsstrahlung cross-section from and Vcut = 10~^ are presented in column 
3. The values in column 4 were calculated with the cross- section from [^], Vcut = 
10^^ and with continuous energy loss due to knock-on electron production. The 
simulation related to column 5 has been carried out under the same assumptions as 
that for column 4 but both knock-on electron and pair production have been treated 
as continuous processes. Finally, column 6 shows the results of the simulation similar 
to that for column 4, but with Vcut = 10~^. The accuracy of the values presented in 
Tables |I| and is of the order of 1%. 

Let us first compare the intensities obtained with different assumptions used in 
the code. The intensities calculated with the muon bremsstrahlung cross-section 
from are systematically higher (by (4-8)%) than those with the cross-section 
from 1^ . The relative difference slightly increases with depth. If the depth-intensity 
relation measured in some experiment is used to obtain the power index of the muon 
energy spectrum at sea level (or the primary spectrum), then, the use of the cross- 
section from 1^ will result in a higher power index (softer muon spectrum) than 
that with cross-section 0]. The difference can be of the order of 0.01. 

The treatment of the knock-on electron production and pair production as con- 
tinuous processes results in the decrease of muon intensities (up to 10% at high 
depths). The relative difference also increases with depth. The power index of the 
primary spectrum which can be evaluated from the underground measurements will 
be lower in this case. The difference in the power index is negligible if only knock-on 
electron production is treated as a continuous process but can be of the order of 
0.02 if both knock-on electron and pair production are not treated stochastically. 

The increase of Vcut from 10"'^ up to 10^^ results in a small decrease of the muon 
intensities ((2 — 3)%). However, the increase of the relative difference with depth 
is small and may not influence much the power index which can be evaluated from 
the measured depth - intensity relation. The increase of power index in this case 
should not exceed 0.01. 

All aforementioned treatments lead to the conclusion that the reasonable cross- 
sections and assumptions used in the code (may be with the exception of the con- 
tinuity of the pair production process) result in the muon intensities which do not 
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Table 1: Vertical muon intensities, I{h), in units cm~'^ s^^ sr"^ for different depths 
in the standard rock: column 1 - depth in km.w.e.; column 2 - bremsstrahlung cross- 
section from [4], Vcut = IQ-^, all processes are stochastic; column 3 - bremsstrahlung 
cross-section from [9], Vcut = 10~^, all processes are stochastic; column 4 - 
bremsstrahlung cross-section from [4], Vcut = 10~^, continuous energy loss due to 
knock-on electron production (e - cont.); column 5 - bremsstrahlung cross-section 
from [4], Vcut = IQ-^, continuous energy loss due to knock-on electron production 
and pair production {e + p - cont.); column 6 - bremsstrahlung cross-section from 
[4], Vcut = 10-2, continuous energy loss due to knock-on electron production. 
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differ much from each other. The difference is comparable with the systematic uncer- 
tainty which can be attributed to the existing experimental data. The difference in 
the power index of primary spectrum which can be evaluated from the experimental 
data using calculated intensities should be less than the errors (including system- 
atical error) of the existing data. However, the stochasticity of pair production is 
quite important for this purpose. 

We would like to note that these considerations are attributed to the single 
muons. The situation can change for multiple muons when the small decrease of 
the muon survival probability can result in the significant decrease of the intensity 
of muon bundles with large multiplicity. In section 4 we present some investigations 
of this point. 

The mean muon energies calculated with our basic assumptions (2nd column 
of Table ^ are higher than those obtained with another cross-section and/or other 
assumptions. The difference is not negligible in the cases of the cross-section from 
1^ and/or the continuity of the pair production process. 

We have carried out also the simulation of muon transport with our basic as- 
sumptions and with Vcut = 3-10-^. We observed no statistically significant difference 
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Table 2: Mean muon energies at vertical, < E{h) >, in GeV for different depths 
in the standard rock (see caption of Table 1). 
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in the muon intensities and mean muon energies comparing with those obtained with 
Vcut = 10-3. We can conclude that our basic assumptions are optimal from the point 
of view of the correctness of the results and CPU time. 

We have checked also the influence of the cross-section difference and Vcut on 
the angular and lateral distributions. The initial sea-level muon spectrum in the 
power-law form with power index 3.7 has been used. With the cross-section from |^ 
we have obtained the mean angular deviation 0.56° and 0.45° for 3 and 10 km.w.e., 
respectively, while with the cross-section from the values are 0.56° and 0.46°. 
Statistical uncertainty is of the order of 0.01°. The mean lateral displacement is 
1.56 m and 1.87 m with the cross-section from for 3 and 10 km.w.e., respectively, 
and 1.59 m and 1.93 m with the cross-section from [^. Statistical uncertainty is 
about 0.01 m. We can conclude that the difference in the cross-sections does not 
influence much the angular and lateral distributions at any depth. A small but 
statistically significant difference of the mean lateral displacements at large depth 
can be attributed to the difference of the mean muon energies at this depth (see 
Table g and Eq. (5)). 

We have carried out also the propagation of muons with Vcut = 

10-2. Actually, 

the increase of Vcut from 10-3 to 10-^ results in an increase of the average step 
between two consecutive interactions by one order of magnitude. Thus, this proce- 
dure allows us to test the dependence of the resulting angular deviation and lateral 
displacement on the step between two consecutive interactions. We didn't find any 
significant difference between the results with Vcut = 10-3 and Vcut = 10-^. This 
proves the correctness of the three-dimensional simulation procedure. 
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4 Muon multiplicity 



We have roughly investigated also the effect of bremsstrahlung cross-section on muon 
multiplicity distribution. We have used the code CORSIKA [|TU| to simulate the 



development of a shower produced by a primary cosmic-ray nucleus in atmosphere. 
We have sampled two sets of 10000 showers produced by protons and iron nuclei of 
lO'* TeV to have a good statistics of muons. 

To transport shower muon component through 3000 m.w.e. of standard rock with 
MUSIC, we have used two cross sections for bremsstrahlung from and [Q. For 
the case of infinite underground detector, we have produced two muon multiplicity 
distributions obtained with different cross-sections. Figures |^ and [1^ show these 
distributions for showers produced by protons and iron nuclei. In Table § the mean 
values < > are presented, their percentage difference <^Maa> <^^kkp> jg ^ ^<y^^ 



Table 3: Mean muon multiplicity at 3 km.w.e. for proton and iron primaries of 
energy 10^ TeV obtained using different bremsstrahlung cross-sections. 



Cross — section 


<N^> p 


<N^> Fe 
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7.89 ± 0.04 


20.96 ±0.07 


i 


1 


8.20 ±0.04 


21.71 ±0.07 



In order to evaluate the possibility that the differences between distributions 
could be caused by statistical fiuctuactions, we have done a test: the probability 
that muon multiplicity distributions obtained with different cross sections come from 
a single distribution is ~ 10%. 

The effect of different bremsstrahlung cross sections on muon multiplicity can be 
underlined using a fit with a negative binomial E^. It is evident, as shown in Fig. 



TT|, that fit to muon multiplicity distribution, obtained with cross-section from 



is practically a translation of the other fit due to different values of < A*"^ >. 

These simple considerations proves that the muon interaction cross-sections can 
influence the muon multiplicity distributions underground and, hence, the primary 
cosmic-ray composition which can be evaluated from the underground measure- 
ments. 

In real experiments, however, the observed multiplicity distribution is driven by 
the primary energy spectrum and the assumed primary composition and a such 
complete simulation to investigate the differences due to the muon cross sections is 
beyond the aim of this paper. We can state only that the effect is quite small (but 
not neghgible) and dominated by other systematic uncertainties. 



5 Conclusions 

We have presented the new three-dimensional Monte-Carlo code MUSIC for muon 
propagation through the rock. The main advantages of this code are: i) the treat- 
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ment of all main processes of muon interaction with matter (bremsstrahlung, in- 
elastic scattering, pair production and knock-on electron production) with the frac- 
tional energy loss v > 10^'^ as stochastic processes; ii) the use of the most accurate 
cross-sections of muon interactions with matter; iii) the accounting of the angular 
deviation of muons due to stochastic processes. We can conclude that the angular 
and lateral distributions of muons at any depth are fully determined by multiple 
Coulomb scattering, while other processes give a contribution which is 1-2 orders of 
magnitudes less. However, the very rare processes of muon scattering to large angles 
(for example, backward scattering) which can be due to inelastic scattering and/or 
bremsstrahlung need more accurate treatment of these processes at large Q^. We 
have studied the dependence of the characteristics of muon flux underground on the 
bremsstrahlung cross-section and assumptions used in the code. We can conclude 
that the treatment of pair production as a stochastic process is quite important 
for the characteristics of muon flux underground. The separation of stochastic and 
continuous energy losses at Vcut = 10^^ is optimal from the point of view of CPU 
time and the accuracy of the results. The accurate knowledge of the cross-sections 
is needed to calculate the intensity of single muons and muon bundles underground. 
Even a small (1-2)% increase of the muon bremsstrahlung cross-section 0] comparing 
with that from results in a significant (6%) decrease of the single muon intensity 
at 3 km.w.e. and in a shift (however, quite small) of muon multiplicity distributions 
toward low multiplicities. This could produce some bias in the interpretation of the 
muon 'depth-intensity relation' and multiplicity distribution measured in an under- 
ground experiment when they are analysed in terms of the primary spectrum and 
composition. 
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Appendix A. Simulation of muon angular deviation 
due to stochastic processes 

To sample the deflection angle of muon due to bremsstrahlung or pair production 
processes we followed the approximations for mean angle < 9'^ >^/^ between the 



incident and scattered muon proposed in [0. 



Let us consider first the muon bremsstrahlung process. Since for both energy 



loss and angular deviation the so-called 'coherent nuclear bremsstrahlung' [|T^ 
is the most important, we have used the parametrisation for < 6'^ >^/^ of [O for 
this process: 
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< 02 >l/2 



^5(1 - v] 



5 '"2/ . , 

mm 



< 0.5 

nax 

>^/2< 0.2,1; > 0.5 



-1/2 ^ ^2 >i/2> 0.2, > 0.5 



(9) 



-0.91 



where = 0.092^-^/3, ^ 0.052^-1^-1/^ ^ 0.22E-°-9^ = 0.26^" , 
n = 0.81i?i/2/(£'i/2 + 1.8), and k^ is related to ^4 by continuity of the function. 
The mean angle between incident and scattered muon due to pair production 



process has been parametrised as 112 



< 02 ^1/2 ^ {2.3 + \nE)E-\l-v)-\v-2me/Efv' 



X 



8.9 ■ 10-^1/^(1 + 1.5 ■ IQ-^E) + 0.032t;/(f + 1), 0.1 



(10) 



111 with- 



where rrie is the electron mass. 

For a given value of v sampled according to the cross-sections from 
out any additional approximation the mean angle < 0^ >i/2 between incident and 
scattered muon has been calculated using the Eqs. (9,10). The spatial angle 9 has 
been chosen from an exponential in the variable with the mean value < 6^ >. 

We note again that we did not use any additional parametrisation for the muon 
inelastic scattering cross-section which has been taken from |T^. For a given value 
of V sampled according to the cross-section integrated over Q2, the muon scattering 
angle 6 has been simulated using the double-differential cross-section. 

For all stochastic processes we have sampled randomly the angle in the plane 
perpendicular to the muon direction. 



The code is available upon the request by e-mail to: antonioli@bo.infn.it or 
kudryavtsev@vaxmw.tower.ras.ru 
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Figure 1: Energy distributions of muons (from a power spectrum with index 7 = 
3.7) at 3 km.w.e. (solid line) and at 10 km.w.e. (dashed line). 
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Figure 2: Angular distributions of muons (from a power spectrum with index 7 = 
3.7) at 3 km.w.e. (solid line) and at 10 km.w.e. (dashed line). 
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Figure 3: Lateral distributions of muons (from a power spectrum with index 7 = 
3.7) at 3 km.w.e. (solid line) and at 10 km.w.e. (dashed line). 
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Figure 4: Comparison of angular distributions obtained with MUSIC (solid line) 
and PROPMU (dashed hne) at 3 km.w.e. 
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Figure 5: Comparison of lateral distributions obtained with MUSIC (solid line) and 
PROPMU (dashed hne) at 3 km.w.e. 
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Figure 6: Angular deviations generated by each mechanism of deflection separately 
at 3 km.w.e. (from a power spectrum with index 7 = 3.7): multiple scattering (con- 
tinuous: Gaussian theory, dashed: Molicrc theory), inelasting scattering (dashed), 
bremsstrahlung (dot-dashed), pair production (dotted). 
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Figure 8: Comparison between survival probabilities obtained with MUSIC (solid 
line) and PROPMU (dashed hne) at 4 km.w.e. and 10 km.w.e. 
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Figure 9: Comparison of muon multiplicity distributions at 3 km.w.e. obtained using 
bremsstrahlung cross sections from [9] (open circles) and [4] (triangles). Muons are 
produced in 10000 proton showers of energy Eo = 10'^ TeV 
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Figure 10: Comparison of muon multiplicity distributions at 3 km.w.e obtained using 
bremsstrahlung cross sections from [9] (open circles) and [4] (triangles). Muons are 
produced in 10000 showers initiated by iron nuclei of energy Eg = 10^ TeV 
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Figure 11: Comparison of the fits to the multiphcity distributions with nega- 
tive binomial functions. Multiphcity distributions have been obtained using muon 
bremsstrahlung cross sections from [9] (dashed hne) and [4] (sohd line) for proton 
and iron primaries. 
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